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Abstract 
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Alternating patterns of small and large amplitude oscillations occur in a wide variety of 
physical, chemical, biological and engineering systems. These mixed-mode oscillations (MMOs) 
are often found in systems with multiple time scales. Previous differential equation modeling 
and analysis of MMOs has mainly focused on local mechanisms to explain the small oscillations. 
Numerical continuation studies reported different MMO patterns based on parameter variation. 
C/5 ' This paper aims at improving the link between local analysis and numerical simulation. Our 

f"*^ . starting point is a numerical study of a singular return map for the Koper model which is a 

prototypical example for MMOs that also relates to local normal form theory. We demonstrate 
that many MMO patterns can be understood geometrically by approximating the singular maps 
with afhne and quadratic maps. Motivated by our numerical analysis we use abstract affine and 
quadratic return map models in combination with two local normal forms that generate small 
oscillations. Using this decomposition approach we can reproduce many classical MMO patterns 
£\J . and effectively decouple bifurcation parameters for local and global parts of the flow. The overall 

^ ' strategy we employ provides an alternative technique for understanding MMOs. 
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Complex oscillatory patterns have been observed in a wide variety of applications. 
Analyzing these patterns from a dynamical perspective has been an active area of 
research for decades. However, several mathematical breakthroughs in the last 15 
years have provided substantial additional insight into phenomena that describe lo- 
^ cal oscillations. In the present paper, we provide a numerical study of the singular 

Poincare map in the Koper model. We demonstrate that many MMO patterns for 
the Koper model can already be understood just using approximations of singular 
limit maps. The results for the Koper model suggest that a local-global numerical 
simulation approach combining normal forms with discrete maps can be effective. We 
show that this abstract approach reproduces many typical MMO patterns that have 
been observed in applications. This methodology aims to close a gap between previous 
numerical studies of MMO patterns and analytical results about local normal forms. 

1 Introduction 

Mixed-mode oscillations (MMOs) are patterns of small and large amplitude oscillations in a time 
series that differ at least by one order in magnitude. They have been observed experimentally in 
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the Belousov-Zhabotinsky reaction in the 1970's and 1980's \39\ [53] and have been encountered 
more recently in a wide variety of different experiments [3T| . I38| . f6T |, I19J . The basic classification has 
been based on counting the number of small oscillations s and large oscillations L so that we can 
symbolically represent an MMO by 

r s j— i t s i t s i+ 1 

■ ■ ■ L j~l L j L j+1 ■ ■ ■ 

where j 6 N is an index. For example, if we have a periodic time series that has 2 large amplitude 
oscillations (LAOs) and then 3 small amplitude oscillations (SAOs) we get . . . 2 3 2 3 2 3 ... or simply 
an MMO of type 2 3 . Systems exhibiting MMOs are often modeled using differential equations 
[HJ [3] . Local bifurcation theory [101 [29] and numerical methods [33 EH] have been developed to 
gain a lot of insights into SAO generating mechanisms [9]. A multiple time scale structure of the 
system is a key component for many local mechanisms. A detailed survey of this theory and its 
applications to particular models has been completed recently [15]. The main findings of many 
numerical studies (see e.g. [351 [731 E21 EZ]) and experiments (see e.g. [391 E3 EQ]) are transition 
sequences of periodic orbits; for example, if we only consider MMOs with patterns of the form 
• • • L S L S ■ ■ ■ such a transition sequence can be represented as follows 

• ■ ■ -> (L pi y*i -> (L P2 y^ -> (L P3 y*s -> . . . (i) 

where p is a control/bifurcation parameter i.e. under variation of a single parameter changing pat- 
terns of MMOs can be observed. To understand patterns of the form ([1]) several approaches have 
been used. The theory of local normal forms has been applied to explain the SAOs and then it 
is usually assumed that the global return mechanism satisfies certain properties (see e.g. [731 110] ) 
so that the local theory becomes applicable or a phenomenological model for the return map is 
proposed \58\ 159] . Another approach is to compute Poincare maps [52] numerically under param- 
eter variation (see e.g. [42 \ [56] 133] ) to explain transitions of MMO patterns or to use numerical 
continuation [46J to subdivide parameter space (see e.g. [33 EH]). These techniques have provided 
tremendous insight into what types of sequences ([!]) can be found in different systems. However, all 
previously mentioned studies vary parameters in such a way that local and global dynamics change 
simultaneously. Here we suggest that to understand which patterns of the form (JT]) occur one also 
has to ask what happens when this parameter coupling is not present. Only in this context one is 
able to distinguish the effects of parameter variation on the local normal form from the variation 
of parameters in the Poincare map. We start by applying this idea in the context of Koper's model 
[43] . For Koper's model the local dynamics is well-understood [15] and SAOs are generated by 
folded nodes (STJ [71] and folded saddle- nodes of type II (or singular Hopf bifurcation, [2S1 115] ) 
which are normal forms for systems fast-slow systems with three variables (see also Appendices 
IA.2IIA.3l for a brief review) . 



Remark: We point out that folded nodes and singular Hopf bifurcation are two possible normal 
forms under the assumptions of fast-slow systems structure and non-degeneracy assumptions for a 
folded critical manifold. Obviously one can also suggest other possible SAO mechanisms [35| 156]. 
However, we have chosen to focus on the Koper model that is well-described locally by the two 
normal forms described above. The main reasons for this choice are that many experimental and 
analytical studies have been found that exhibit folded nodes and/or singular Hopf bifurcation (see 
the review [15] for a list systems with folded nodes and singular Hopf bifurcation) . Furthermore, it 
has recently been shown that both mechanisms also relate to delayed Hopf bifurcation [50] which 
has been proposed as another SAO mechanism. 
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The global return mechanism for MMOs in the Koper model is provided by a cubic relaxation- 
oscillation mechanism [44, 68J that has already been investigated by van der Pol in the 1920s |13U14|. 
Here we provide numerical computations of the global Poincare return map as a composition of 
several maps in the singular limit of perfect time scale separation. These calculations reveal that the 
return map can be surprisingly regular. Using affine and quadratic approximations to the singular 
maps we investigate MMO patterns and find that the approximations suffice to understand MMO 
sequences observed in extensive numerical continuation. Motivated by these results we combine 
two local normal form ODEs with abstract linear and quadratic maps to study MMOs. It is shown 
that classical sequences of the form ([1]) as well as chaotic MMOs can be easily generated in this 
framework. In particular, it is easy to design MMO patterns and to understand the differences in 
local and global parameter effects. We point out that this study also contributes to closing the gap 
between numerical simulation and local normal forms by reproducing several of the MMO transi- 
tion sequences observed by a simultaneous local and global parameter variation in the Koper model. 

The paper is structured as follows. Appendix [A] contains the necessary background for readers 
not familiar with fast-slow system and MMO generating mechanisms in these systems. The main 
part of this paper starts in Section [2] where the Koper model is introduced and its basic properties 
are reviewed. In Section [3] the global singular return map for the Koper model is decomposed into 
several more tractable flow maps using numerical simulations. In Section [4] the maps are approxi- 
mated by affine and quadratic map models; Appendix[B] contains a discussion of the approximation 
error. In Section [5] the global aspects of MMOs in the Koper model are analyzed using the flow map 
models. In Section [6] we consider a standard local-global decomposition of the MMO generating 
mechanisms. The key point is that we suggest to separate the parameter dependencies for the local 
and global models. We combine a global return map model with local SAOs induced by folded 
node and singular Hopf normal forms. We conclude with a brief outlook, describing the wider 
applicability of our approach, in Section [7j 

2 The Koper Model 

One version of the Koper model for MMOs is given by 

e\x = y — x 3 + 3x, 
y = kx - 2(y + A) + z, (2) 
z = e 2 {X + y-z), 

where (k, A) are the main bifurcation parameters and (61,62) are the singular perturbation param- 
eters. The equations were first studied as a two-dimensional model by Boissonade and De Kepper 
[7J modeling a prototypical chemical reaction. Koper [45J added a third variable to a planar system 
and used numerical continuation techniques [2S1 El] to study MMOs [45 1 . It is very important to 
note that equations similar or equivalent to (J2|) have been proposed many times independently by 
several different research groups [251 ESI HS1 Ell S3 EH] as a "canonical" , "minimal" or "typical" 
model for MMOs. The version ([2]) of Koper's model was proposed by the author and co-workers in 
|15j : it is obtained by a coordinate transformation of Koper's original model and has the symmetry 

(x, y, z, A, k) H> (-x, -y - z, -A, k) 

which allows us to restrict to parameter regions with A > or A < without loss of generality. Here 
we shall only review the local bifurcation structure briefly an introduce the necessary notation; a 
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detailed local fast-slow systems analysis of ([2]) can be found in |15| . We also point out that the 
terminology reviewed in Appendix [A] will be assumed from now on. 

If < 61,2 *C 1 holds then ([2]) is a three time-scale system. We shall focus on the case £2 = 1 
and < ei =: e <C 1 in which case we have one fast variable x and two slow variables (y, z). The 
critical manifold is 

C = {(x, y,z) £R 3 :y = x 3 - 3x =: c(x)}. 
The typical cubic (or S-shaped) structure splits the critical manifold into several parts 

c = c a >- U F_ U C r U F + U C a ' + 

where C a ~ := C n {x < — 1}, C a ' + := C D {x > 1} are normally hyperbolic attracting, C T : = 
Cn{ — l<x<l}is normally hyperbolic repelling and 

F_ : =Cn{x = -l} = {(-l,2,z)} and F + := C n {x = 1} = {(1, -2, z)} 

are fold curves curves of the critical manifold. MMOs can easily be observed in simulations; see 
Figure CD The desingularized slow subsystem is 

x = kx — 2{c(x) + A) + z, 

z = (3 x 2 _3)(A + c(x)-z). (3) 




Figure 1: The parameter values for the simulation are (e, k, A) = (0.01,-10,-7). The time series 
for a l x l 2 MMO for the variable x is shown on the left and the phase space trajectory is shown on 
the right; the critical manifold Co is shown in grey. 

There are two folded singularities 

p± = (±l,2A=F(4 + fc)). 

By symmetry we shall only focus on the folded singularity p+. The linearization of the desingular- 
ized slow flow p + is 

i)-U%s)(sws) 
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where we already set k = — 10 which will fixed from now on. Note that in this case {X = 0} 
corresponds to the fold line F + and p+ is located at the origin. The eigenvalues of A + are 

a w (X) = -5 + V-23 - 6A and a s (X) = -5 - V-23 - 6A, 

with associated eigenvectors 

E W (A) = ^ 5+V-23-6A J an d E S (A) = f 5-V-23-6A j . 

Therefore p+ is a folded saddle for A < —8, a folded saddle-node of type II (FSN II [67l [I5j) for 
Afsn ii = —8 and a folded node for A E (—8,-23/6). At A n j = —23/6 the transition from a folded 
node to a folded focus occurs. The singular Hopf bifurcation for the full system occurs 0(e) away 
from Afsn ii- It is supercritical and the stable global equilibrium q loses stability at this point. 
Therefore the interesting parameter region for MMOs is 

A 6 (A F sn ii, Kf) = (-8, -23/6), k = -10. 

Note that this parameter region represents a typical one-parameter MMO sequence |45[ [To] . The 
important eigenvector for global returns is S s associated to the strong primary canard 7 S as it 
bounds the rotational sectors lying on C a,+ . The x-component of S s lies, for the scaling we 
have chosen, between Sf(— 8) = oo and £J(— 23/6) = i. Therefore the rotational sectors [10] 
that subdivide the funnel region are given by a convex cone with opening angle between ^ and 
cos" 1 (25/26). 



3 Return Maps - Decomposition 

Our goal is to analyze the structure of the global singular return map. Instead of using the standard 
approach of computing the Poincare map between two fixed sections |64l I3"T| I5"2] we are going to 
decompose the map according to fast-slow systems theory (see e.g. [361 EH EE] for this approach). 
We fix e = and recall that k = —10. Then we focus on A as the primary bifurcation parameter. 
In this case the folded node p + and the unique equilibrium q account for the SAOs. Observe that 
global returns to a neighborhood of p + can be decomposed. See Figure [2] for an illustration of the 
one-dimensional singular maps we are going to define: 

(a) Trajectories can reach the fold line F + at a jump point and follow the fast flow to the drop 
curve L a ~ := CC\ {x = —2}. Then trajectories follow the slow flow induced by ([3]) to F_ and 
jump to the drop curve L a,+ := C fl {x = 2}. We denote this map by 

mj : F + -> L a ~ -)• F_ -»• L a ' + 

where j indicates that we consider a regular jump. We denote the intermediate map by 
m a _ : L a '~ —J- Observe that if we parametrize the domain and range by z then the 

intermediate map m a - is the only non-trivial component of the map rrij and the other parts 
of mj are the identity with respect to z. 

(b) Trajectories can flow into the folded node p+. Suppose we consider trajectories tracking the 
part of the strong canard 7 S contained in C r . These trajectories jump at some point from 7 S 
to C a ~ and flow into F_ before jumping to L a ' + . Denote this map by 
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where / indicates a jump forward (or away) singular canard orbit; again observe that only 
the part C a ~ — > F_ is non-trivial with respect to z. 

(c) Trajectories tracking the strong canard "y s C C r can also jump at some point from 7 S to C a,+ 
and flow into F + . It will be advantageous to terminate this map at a line := Cf){x = 1 + 
for some fx > sufficiently small. Then we have a map 

where b indicates a jump backwards (or back) singular canard orbit. 

(d) There is also a map induced by the slow flow on C a ' + starting from the drop curve L a ' + 
towards the fold line 

m a ,+ : L a ' + -> L M 

(e) The linearization ()4]) at the folded singularity p + can be used to define a flow map in the fold 
region 




Figure 2: Illustration of the singular map decomposition; parameter values are (e, k, A) = 
(0, —10, —7). Definitions of all maps and domains are given at the beginning of Section [3l Here 
we show: the critical manifold C (grey), the strong canard 7 S C C m (green) and its projections to 
C a,± (dashed green), the fold lines F± (black) and their projections L a,=F (dashed black, fx = 0.1), 
the line (yellow) and the folded node p+ (blue circle). Examples for the maps rrij (red, regular 
jump), m a ^ (blue, flow towards F + ), nif (magenta, jump forward canard) and nib (magenta, jump 
backward canard) are displayed as well. 

Figure [3] shows representatives of the maps rrij, m aj+ , and rrif for A = —7 with respect 
to the variable z and also the associated slow flows. The main observation is that the maps are 
surprisingly regular. 
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Figure 3: Singular maps for (k, A) = (—10, —7) with respect to the z- variable i.e. the horizontal axis 
shows z = Zi n and the vertical axis shows z out = rriK(zin) for K £ {j, (a, +), b, /}. The insets (grey 
background) illustrate the phase space flow on the attracting critical manifolds C a,± associated to 
the maps mx] we show only every tenth trajectory in the computation of mx- The color coding is 
the same as in Figured 

4 Return Maps - Modeling 

In this section we are going to discuss the modeling of the maps computed in Figure O a discussion 
of the approximation error as well as the error for e > is given in Appendix [Bj Figure [3] motivates 
considering afhne and/or quadratic maps. The map rrij seems to be close to an afhne map which 
is due to the very simple regular slow flow from L a ~ to F_ ; see also Figure (2) Similarly, we 
propose to model the map m aj+ by a quadratic map. The maps induced from the projections of 
the strong canard 7 S C C r onto C a ^ are multi-valued when parametrized with respect to z due to 
the fold structure of 7 S ; see Figure [2j With another parametrization we expect that and m f 
are generically single-valued by uniqueness of solutions for the desingularized slow subsystem. The 
parametrization with respect to z is very convenient. We propose to make the following ansatz: 

mxiz) = c 2 (A)z 2 + ci(X)z + co(A) 

for each map uik with K G {j, (a, +), b, /} where the coefficients Co,i,2(A) are to be determined. 
We are going to illustrate the procedure for finding the coefficients for mf and just state the results 
we obtained for the other three maps. The ansatz is that mj can be decomposed as follows: 

r c{ u (\)z + 4 u (x) if zi u m (\) < z < zLux), 

m f (z) = < 4 l (X)z 2 + c{\X)z + c f l (X) ifi(A) < z < z£L x (X), (5) 
I, undefined otherwise, 

where we impose continuity at the shared boundary point rrif(z-^ in ) = rnf(z n 1 lin ). See Figured] for 
an example. In Figure H] the upper part of mj is approximated by an affine map and the lower part 
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by a quadratic. The computation of the approximation error in Appendix [B] for A G (Afsn ii, A n /) 
shows that for each fixed value of A the affine and quadratic models provide an approximation on 
the order of 10 -2 of the singular maps obtained via numerical integration of slow flow trajectories 
on a fine mesh. 




Z min Z min 



Figure 4: Singular map 77i f with. z ou i — ^/(^in 

). The computed map is shown as a dashed grey 
curve and the approximations are shown in solid black (affine for upper part and quadratic for lower 
part). The bounds of the domains for each part of the map are marked as well (dotted vertical 
lines). 

As a next step we investigate all functions depending on A in ([5]). The boundary Zmax is given 
by the folded singularity p+ so that Zmax(X) = 2A + 6. We also know from the definition of © that 



Jl 



The other functions of A can only be approximated numerically due to the nonlinear 
slow flows on C r , which defines j s , and on C a ~ , which defines the map to -F_. Figure [5] shows 
numerical computations of the unknown functions of A in the definition of mj in ([5]) . 

Several observations can be made from Figure [5] and the previous remarks. All the domain 



boundaries z^ in (X), Zmax{X), z minW an d Zmax(X) are very regular and seem to depend almost 
linearly on A. The coefficients of the linear and quadratic polynomials have substantial nonlinear 
dependencies on A for the entire range A E (Afsn ii, A n /)- This implies that although affine and 
quadratic maps can be very good approximations at fixed parameter values it will be more difficult 
to analyze the global return maps inducing MMOs as parameter-dependent families. For the other 
maps m&, rrij and m a ,+ we propose the following approximations: 



J u 



Jl 



Jl 



m b (z) 
rrij(z) 

m a,+ (z) 



c^(X)z 2 +c\ u (\)z + 4 u (\) 
c\ l (\)z + c b l {\) 
undefined 

c5(A)z 2 + C ?(A)z + cS(A), 



^ z minW 

Ji >(A) 



< z < z max (X), 
w (A), 



if z min(^) — z — z ma • 



otherwise, 



(6) 

(7) 

(8) 
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-8 -6 -4 -8 -6 -4 -8 -6 -4 




-8 -6 -4 -8 -6 -4 -8 -6 -4 



Figure 5: Horizontal axes are A and vertical axes are the respective coefficients e.g. top left figure 
shows z^ in (\). The only relevant part for the definition of ([5]) that is not shown is Zmax which is as 
regular (almost linear) as the other parts of the domain boundaries for rrif. The dots are computed 
points and the curves provide polynomial fits (quadratic=blue, cubic=green and quartic=red). 



where we impose continuity at the shared boundary point for mj i.e. mb( z min) = m b( z min)- As 
a next step we are going to calculate the map for the linearized desingularized slow flow near p + . 
The intersection of the eigendirection of T, s with = {x = 1 + /i} is easily calculated as 

(1, 2A + 6) T + L, = (1 + fi, 2A + 6 + //(5 - V-23-6A)) =: (1 + », z"(A)). 

Hence all trajectories that arrive at with z > z^{\) will stay in the funnel and reach p + while 
trajectories for z < z^(X) will first reach the fold line F- and jump to L a ~. To see where on F + 
the last class of trajectories ends up we could just solve Note however that there exists an 
approximation for z < z^(X) that just amounts to projecting (/x, ^(0)) parallel to S s onto {X = 0} 
which is given by 

Ou,Z(0))h+ (o,z(o)- £ 
Therefore we get the local representation for the map m s in (X, Z)-coordinates 

m loc (z) _{ 2A + 6 if^>^( A)) 



If z is the coordinate obtained in original coordinates without linearization then 

m s (z) = 



2A + 6 ifz>^(A) 
z-$e if«<«"(A) 

where the error is 0(/i) as \i — > 0. With the different maps available we can proceed to analyze 
how they can be used to explain the global returns that generate LAOs. 
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5 Mixed-Mode Oscillations 



Throughout this section we work with the polynomial approximations to the maps mi\ that have 
been derived in the last section. The first question we shall consider is what happens to trajectories 
that do not follow the canard 7 S C C r when arriving at p+ or which land outside of the funnel 
region. The relevant map for this purpose is 

(m a ,+ o rrij) : F + n {z < 2A + 6} -> V* (9) 

We are interested when part of the domain of Q is returned inside the funnel so that {m a ^ + o 
rrij)(z) > z fi (X). Figure M shows the map ([9]) for three different values of A. We observe that 
closer to the folded saddle-node of type II (i.e. near the singular Hopf bifurcation) trajectories that 
arrive outside the funnel on F + can get mapped back into the funnel under ([9]). For A = —6.5 
in Figure [6] we observe that no trajectories can return into the funnel and that the return map 
(m s o m aj+ o rrij) will have a stable fixed point since [i is small and hence the projection m s will 
preserve the intersection with the diagonal. 




-11 -10 -9 -10 -9 -8 -9 -8 -7 



Figure 6: Map {m ai+ omj)(z), approximated by (|7|) and ([8]) with k = —10 and /i = 0.1. Horizontal 
axes are input z-coordinates on a domain z G ((2A + 6) — 2, 2A + 6) C F + and vertical axes are 
o rrij)(z) (think black curves). The location of the folded node funnel region z^(X) is shown 
by horizontal dashed black lines and the diagonal is indicated by the thick grey line. 

Hence we can consider several quantitative questions: 

1. For what values of A do trajectories from outside the funnel re-enter it? 

2. When does the map (m s o m a) + o rrij) have fixed points? When does the fixed point coincide 
with the folded node p+? 

3. How are trajectories mapped into the funnel? More precisely, what is the dependence of the 
distance 5 to the strong singular canard 7 S n C a ' + upon varying A? 

A trajectory starting for z < 2A + 6 will re-enter the funnel after one global return if and only 

if 

{m a& o mj ){z) = c^(A)(c{(A)z + ^ (A)) 2 + c?(A)(c{(A)z + C J (A))+cS(A) 

= c a 2 {c\) 2 z 2 + (2c§cj + cfcj) z + (4) 2 c a 2 + cl4 + eg < *"(A) 

By monotonicity of ([9]) on the required interval (see Figure [6]) we can just pick the folded node 
z = 2A + 6 and determine when the condition fails; this yields the critical parameter value at which 
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not all trajectories near p + return to the funnel in one iteration. We find that the parameter value 
at which p + gets returned to the boundary of the funnel is A = A r ~ —6.7887. Next, we consider 
the fixed points of (m s o m a ^ + o rrij). Those points correspond to candidates representing relaxation 
oscillations. We find that at A = A r a stable fixed point appears for the map (m s o m aj+ o rrij). 
Therefore we find that a transition to relaxation oscillations occurs near A r for the full system and 
e sufficiently small; this can be confirmed by numerical continuation |15| . Note that the bifurcation 
that creates the fixed point occurs at the boundary of the domain of (m s o m a; + o rrij). 

As a next step we consider candidates that follow the canard 7 S n C r i.e. we consider the 
maps rrif and m&. We start with m& which represents medium-size canard- induced oscillations if 
trajectories from the domain of m& re-enter the funnel after one iteration step. Figure [7] plots three 
examples of the map m&. The closer the parameter values are to the folded saddle- node of type II 
at A = —8 the larger is the part of 7 S n C r that returns inside the funnel. The closer we are to 
relaxation oscillation at A = A r the more of j s n C r gets mapped outside the funnel. 




-9.5 -8.5 -8.5 -8 -7.5 -7.5 -7 -6.5 



z z z 

Figure 7: Singular maps with respect to the z-variable i.e. the horizontal and 
the vertical axis shows mj, = m^Zm). The horizontal dashed line indicates the funnel boundary 
z^(A); here fi = 0.1. (a) A = -7.5, (b) A = -7 and (c) A = -6.5. 

Note that near A = — 8 with A > — 8 we must always have some part of 7 S n C r near p + that 
does get mapped outside the funnel since the opening cone angle of the funnel region is less than 
5; see Section [2j Therefore there is always one part inside and one part outside the funnel for jump 
back canard orbits. Orbits in the full system that follow 7! for an 0(l)-time on the slow time scale 
and get mapped back to C a ' + via perturbation of m& represent intermediate oscillations. 

For the map m/ we immediately consider m a) + omj to see how jump forward canards get 
returned relative to the funnel. Figure [5] shows that there is a very rapid transition from jump 
forward canards that end all in the funnel for A = —7 (Figure EUa)), a splitting of jump forward 
canards with respect to the funnel (Figure EJb)) and all jump forward canards outside the funnel 
for A = —6.6 (see Figure[8^c)). Let us consider the case when the entire jump forward canards end 
up in the funnel. This can be interpreted as a global MMO generating mechanism via canards. 
More precisely, a trajectory of the full system can make small oscillations near a folded node, follow 
C r n 7 S closely producing an intermediate oscillation and then return into the funnel. This provides 
a mechanism to transition small loops into large ones via canards. The closer we get to A = A r the 
more excursions outside the funnel occur which means that in this region we expect more mixed 
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Figure 8: Singular maps with respect to the z-variable i.e. the horizontal axis shows z — Zi n and 
the vertical axis shows m at+ o mf = (m at + o mf)(zi n ). The horizontal dashed line indicates the 
funnel boundary z M (A); here /x = 0.1. (a) A = —7, (b) A = —6.8 and (c) A = —6.6. 



behavior of MMOs of type L s with L > 1. It is also expected that period-doubling bifurcations of 
the return map can explain transitions between regions of different LAOs. Since resonances for the 
eigenvalues [71] of the folded node are fewer near A r we also expect s to decrease if we increase A. 
Hence we find that MMO sequences near a singular Hopf bifurcation will produce patterns with 
s> 1 and small L while away from the singular Hopf L s patterns with L ~ s are more likely to 
occur. All these findings agree with numerical continuation results in [45(. I15j. 

Therefore one main conclusion from the numerical simulations considered here is that the sin- 
gular limit decomposition is already sufficient to explain many MMO transition sequences. Indeed, 
in the singular limit we could already identify the local normal forms (see Section [2] and [15]) and 
here we calculated a decomposition of the global return map. The main point is that we have used 
a different, and easily implementable, numerical technique to understand geometrically many of 
the MMO patterns that have been found using extensive numerical continuation runs |45j . 

6 A Local-Global Model 

We have seen that the global singular return maps for the Koper model are very regular and can 
often be described as affine or quadratic maps. The only feature of the global returns that is 
complicated to describe are canard orbits that follow the strong canard 7 S n C r . These orbits 
describe intermediate oscillations i.e. orbits that, under parameter variation will grow to a large 
relaxation loop or decay to a small oscillation. However, many MMO transitions can be understood 
without these orbits as shown in the previous section. Hence it is natural to ask what happens if 
we do not consider these intermediate orbits and look at a simulation model for MMOs containing 
local and global maps. The local description of this model is chosen as a flow map for a folded node 
or a folded-saddle node ODE normal form; see Appendix [Al We assume without loss of generality 
that the folded singularity is located at the origin (x,y,z) = (0,0,0). For the local dynamics we 
use the normal forms (|15p and (|18p . Recall that the critical manifold of both normal forms is 

C = {(x,y,z) GR 3 :y = x 2 } 

It is attracting for x > and repelling for x < and we denote the two branches of Cq by Cq and 
Cq. The associated attracting slow manifold provided by Fenichel Theory is 

C? = {(x,y,z)eR 3 :x = h a e (y,z)} 
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where the map h® is given by the implicit function theorem and h${y, z) = y/y. Define two sections 



£1 

S 2 



{(x,y,z) G 
{(x,y,z) e 



-k 2 } 



for suitable fixed kj > 0, kj = 0(y/e) with j = 1,2. The choice of scaling 0{y/e) is prescribed by 
the fact that outside of a neighborhood of size 0{y/e) of the origin Fenichel Theory applies. Define 
a map 

m 12 : Ei -> S 2 (10) 
by the flow map of (j!5j) or f)18|) . Note that the sections Sj are naturally parametrized by the 



coordinates (y,z). The global return map m 2 i : X 2 

m 2 \{y,z) = 



Ei will be modeled as follows: 



k\ 
m{z) 



+ € 



an 
a-21 



«12 
022 



+ 



b-2 



+ 0(e 2 ) 



(klm(z)) T + e[A(y,zf + b} + 0(e 2 



(11) 



where m{z) = m 2 z 2 + m,\z + mo and we require that the matrix A is invertible. Note that we can 
make several further choices e.g. we could decide to include higher-order terms or to assume that 
m{z) is modeled as an affine map and set m 2 = 0. The map m 2 i has to satisfy a further constraint 
if we assume that all trajectories approach the origin exponentially close to the slow manifold C"; 
this requires 

fc 2 Ve = {hi o m 2 i)(y,z). 

Another constraint to generate MMOs is that the global map to 2 i maps some part of its domain 
close to the perturbation of the folded node funnel region. Although the model (llip has formally 
nine free parameters TOj, ajk, b\ we can also view m 2 i as an 0(e)-perturbation of the leading order 
term which has only two or three parameters depending on the choice of model (m 2 = or m 2 ^ 0). 
Hence the description is low-dimensional, explicit and decouples the local and global bifurcation 
structure of the problem. To illustrate the effect of global bifurcation parameters we numerically 
investigate two typical MMO sequences, one for each local normal form with fixed local parameters. 

Remark: In the following, we are going to visualize the LAOs in a time series for the dynam- 
ical system defined by mi 2 and m 2 i by inserting a large amplitude oscillation at fixed amplitude 
whenever the map m 2 i is applied. 

For the folded node (|15h we fix the parameters (e, fi) = (0.01, 0.006). Since 2/c + l<^<2A;-|-3 
for k = 82 we know from folded node theory that there will be k + 2 = 84 canards [10\ 115] . The 
theory also predicts that the maximum number of small oscillations is k + 1. In Figure [9] we varied 
the global return mechanism to demonstrate that we can systematically reach sectors near a folded 
node with a sub-maximal number of oscillations and different MMO signatures. The global return 
map is chosen as the lowest order approximating linear map with m{z) = O.lz + mo and A = 0, 
6 = 0. The parameter mo is viewed as the main bifurcation parameter and controls the entry of 
trajectories to the different folded node rotation sectors. We find the following MMO signatures: 
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Maximal MMO signatures can also be obtained but many of the small oscillations will be at an 
exponentially small scale due to the contraction towards the weak canard [15]. Observe that we 
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Figure 9: MMOs generated by the dynamical system with local dynamics (I15p (parameters e = 0.01 
and [i = 0.006) and global dynamics m(z) = O.lz + m,Q with k\ = \fi. The transition from 
local to global dynamics has been applied when x < —\fe. Trajectories have been started at 
(x,y,z) = (^/e,e,0.15). 

have efficiently de-coupled the local parameter dynamics from the global parameter dynamics. 

The second simulation focuses on the singular Hopf normal form (|18j) . We fix the parameters 
(e, v, a, b, c) = (0.01, 0.01, 0.5, —1, 1). For the global return map we again consider the singular limit 
of 77i2i with parameter 7712 = 0, mi = 0.1 and primary bifurcation parameter mo- It is easy to check 
that the equilibrium near the fold curve is located at q = (x eq , y eq , z eq ) w (— 6. 63729 x 10~ 3 , 4.40537x 
10" 5 , -6.63729 x 10" 3 ). The equilibrium is a saddle-focus with one-dimensional stable and two- 
dimensional unstable manifold. We are in the regime where the SAOs are generated/amplified via 
the singular Hopf mechanism. Figure [T0l shows the typical SAOs with increasing amplitude as they 
approach W u (q). We find the following MMO signatures: 
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m = 0.01 




Figure 10: MMOs generated by the dynamical system with local dynamics (I18p (parameters 
(e,u,a,b,c) = (0.01,0.015,0.5,-1,1)) and global dynamics m(z) = O.lz + mo with k\ = yfe. 
The transition from local to global dynamics has been applied when x < —\fe. Trajectories have 
been started at (x,y,z) = {\fe, e, 0.15). 
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It is important to note that the number of SAOs for the singular Hopf bifurcation is not only 
influenced by the folded node but also by the relative positions of the invariant manifolds of q [151 . 
In particular, Guckenheimer [29j points out that the one-dimensional stable manifold W s (q) seems 
to interact in an intricate way with MMO trajectories. Our decomposition approach is well-suited 
to investigate this dependency further once the local unfolding of the singular Hopf bifurcation 
is better understood [32J . Also for the singular Hopf bifurcation we have been able to reproduce 
a typical MMO sequence without varying the local parameters. Extensive additional numerical 
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simulation showed that it is difficult to produce periodic sequences of MMOs of the forms 

L s , with L > 1, s > 1 and L^L^ 2 ■ ■ ■ (12) 

by varying further parameters in the map m(z). These simulations confirm parts of the incomplete 
theory for MMOs in three dimensions |47} HH] which predict the limited number of MMO patterns 
for three time scale systems. Therefore we conjecture that higher-dimensional return maps are 
more likely to account for more complicated MMOs of the form (I12|) . 




Figure 11: Coordinates (x,y,z) before the global map is applied; the horizontal axis shows 
the return number i.e. I s * application of the global map, 2 nd application, etc. The entire or- 
bits are generated by the dynamical system with local dynamics (|18|) (parameters (e, v, a, b, c) = 
(0.01, 0.015, 0.5, —1, 1)) and global dynamics 771(2;) = 3z 2 + 0.2z — 0.8mo with k\ = y/e. The transi- 
tion from local to global dynamics has been applied when x < —y/e. The times series of the returns 
shows typical chaotic non-periodic behavior. 
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T 

Figure 12: Subset of the time series in (x, r) variables associated to the returns in Figure [TTJ 
Irregular oscillations are observed with 1° and l 1 components. 

Chaotic MMO signatures can be produced easily using a suitable quadratic map with 7712 ^ 0. 
Figures [TT] and [12] illustrate an orbit obtained from the dynamical system of the singular Hopf 
bifurcation with global returns generated by the map m(z) = 3z 2 + 0.2z — 0.8; the irregular be- 
havior of the global returns in Figure [TT1 suggests that this orbit is chaotic. It is well-known that 
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systems with two slow variables and one fast variable with S-shaped critical manifold can be chaotic 
[441 l36j l34~t 157] . Koper [32] observed chaotic regions in parameter space in his original analysis of 
([2]); transitions of MMOs to chaotic sequences can also be observed in many other models [15] . 
As shown above, our model is also able to reproduces this aspect of typical MMO models. We 
conclude that our modeling approach reproduces the main dynamical features and decouples the 
global parameter dynamics from the local parameter dynamics. 



7 Brief Outlook 

The strategy and methods we presented in this paper apply, in principle, to any system where the 
MMO mechanism can be decomposed into a local part that generates the SAOs and a global return 
map. For folded nodes and singular Hopf generated SAOs, the overall dimension can be arbitrary. 
Indeed, it has recently been shown by Wechselberger [72] that the local theory in Appendices IA.21 
IA.3I extends to systems with m > 1 fast and n > 2 slow variables. The technique is to use a center 
manifold reduction to get into the situation (m,n) = (1,2). Then we can still compute singular 
return maps as we still have the three important one-dimensional curves that are analogous to L^ 1 , 
7 S and F + in the Koper model. We have resolved the map for the Koper model in more detail using 
the drop curves L a=t . However, we could just compute m aj + o rrij and m a:+ omj- as single maps for 
another system or adapt the finer global decomposition to the fast-slow geometry of the problem. 
Moreover, it is very important to point out that a center manifold reduction has been already used 
in a four-dimensional system with MMOs generated by folded nodes and singular Hopf bifurcations 
|llj . Our methods apply verbatim to the resulting three-dimensional system obtained in |llj . 

It is expected that the return maps for other systems can be more complicated. For example, 
just consider the situation for the Koper model but insert several non-trivial slow subsystem hy- 
perbolic attractors on C a_ . Then the maps rrij and mt may even have gaps since orbits can get 
trapped on persisting attractors on C" _ . Computing singular maps for several well-known MMO 
models [15] and analyzing their structure is an interesting project but is beyond the scope of this 
paper. 

Acknowledgment: I would like to thank two anonymous referees for valuable comments that 
helped to improve the focus and exposition of the paper. 



A Background Review 
A.l Fast-Slow Systems 

We are only going to recall the basic definitions and results about fast-slow systems. There are 
several standard references that detail many parts of the theory [301 [3D EH 1331 E7J H5j [H E]. A 
fast-slow system of ordinary differential equations (ODEs) is given by: 



ex = ep = f(x,y), 

v = = 9(x,y), 



% (13) 



where x 6 K m are fast variables, y £ M n are slow variables and < e <C 1 is a small parameter 
representing the ratio of time scales. Equation (|13f) can be re-written by changing from the slow 
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time scale r to the fast time scale t = r/e 

x = % = f( x > y)i (ia) 

y> = % = eg(x,y). { ) 

The singular limit e — > of (I14p yields the fast subsystem ODEs parametrized by the slow variables 
y. Setting e — > in (fT3|) gives a differential-algebraic equation (DAE), called the slow subsystem, 
on the critical manifold C := {f(x, y) = 0}. Concatenations of fast and slow subsystem trajectories 
are called candidates. 

A subset S C C is called normally hyperbolic if the m x m total derivative matrix (D x f)(p) is 
hyperbolic. A normally hyperbolic subset S is attracting if all eigenvalues of (D x f)(p) have negative 
real parts for p 6 S; similarly S is called repelling if all eigenvalues have positive real parts. On 
normally hyperbolic parts of C the implicit function theorem applies to f(x,y) = providing a 
map h(y) = x so that C can be expressed (locally) as a graph. Fenichel's Theorem pU [40J, [691 [74] 
states that a compact normally hyperbolic submanifold S = Sq of the critical manifold C perturbs 
for e > sufficiently small, including stability and flow properties, to a slow manifold S e . 

A trajectory is called a maximal canard if it lies in the intersection of an attracting and a 
repelling slow manifold. Canards were first investigated by a group of French mathematicians 
[SJEUHUE] usm S nonstandard analysis. Later also asymptotic [231 El [M] an d geometric [22" 1 l4"9 l [67] 
methods have been developed to understand canard orbits. 



A. 2 Folded Nodes 

Normal hyperbolicity can fail in several ways. Here we briefly review the basic properties of two 
such situations [15]. A non-degenerate fold point p E C is defined by requiring that f(p) = and 
{D x f){p) has rank m — 1 with left and right null vectors w and v so that w ■ [(D xx f)(p)(v,v)] ^ 
and w ■ [(D y f)(p)] ^ 0. The set of fold points forms a manifold of codimension one in the Tri- 
dimensional critical manifold C. If m = 1 and n = 2 the fold points generically form a smooth 
curve that separates attracting and repelling sheets of the two-dimensional critical manifold C. 

Two standard generating mechanisms for small oscillations near fold curves of the critical man- 
ifold will be considered in a normal form setup. Br0ns, Krupa and Wechselberger |67J [10] consider 
a normal form 

ex = y — x 2 , 

y = -(n + l)x - z, (15) 



z 



2 - 



where x is the fast variable, (y, z) are the slow variables and n is a parameter. The critical manifold 
for ([T5j) is C = {y = x 2 } with a line of fold points F = {x = 0,y = 0}. F decomposes the critical 
manifold C = C r U F U C a where C r = CD {x < 0} is repelling and C a = CD {x > 0} is attracting. 
Differentiating y = x 2 implicitly with respect to r gives y = 2xx. Therefore the slow flow is 



x 

z 



-{p+l)x—z 
i 

2 ' 



„ 2. . (16) 



Rescaling time by r i— )■ 2xt reverses the direction of the flow on C r and yields the desingularized 
slow flow 

x \ _ ( — (/i + 1) — 1 \ / x 



Z J V /* on*' (17) 

:A 
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The desingularized slow flow has an equilibrium point at the origin = (0, 0) £ F called a folded 
singularity. The eigenvalues (A s , \ w ) = (—1, —fi) of Aq determine the type of the folded singularity. 
It is a folded saddle for fx < 0, a folded node for /x > and a folded saddle-node of type II for 
H = JBTJ [15] . We restrict to the folded node case and \x E (0, 1) here. Then A s is associated to 
the strong eigendirection 7 Si o and X w is associated to the weak eigendirection 7^0- The extension 
of 7s, (7tu,o) under the slow flow is referred to as the strong (weak) singular canard. Trajectories 
in the funnel region bounded by 7 Sj o and F can pass from C a to C r ; see also [67J [10] . 

The singular canards 70, s and 70 w perturb to maximal canards 7e,s and r y e w that lie in the inter- 
section of the two slow manifolds C° D CJ [67] . If N then there are further maximal canards 
arising as intersections of Cf Pi CJ, called secondary canards [71]. In particular, the attracting and 
repelling invariant manifolds twist around each other [301 128] . The number of twists of a trajectory 
in the fold region can be predicted using its distance 5 relative to the strong singular canard and 
by the value of fi [TO]. We agree to the convention that 5 > indicates a trajectory entering the 
funnel region, 5 = describes the strong canard and for 5 < we are outside of the funnel. The 
twists can cause the SAOs of an MMO. 

A. 3 Singular Hopf 

Note carefully that the normal form (|15f) has no global equilibrium point for \i £ (0, 1). However, 
in many applications a Hopf bifurcation occurs near the onset of MMOs [15] which suggests to 
consider the possibility of a global equilibrium point passing through the folded node region. In 
particular, one has to add higher-order terms to the equation for z in (|15p . Augmenting these terms 
it is well-known that the global equilibrium can undergo a Hopf bifurcation at an 0(e)-distance 
from the fold curve. This scenario is also been referred to as singular Hopf bifurcation [8j [29] since 
the pair of complex conjugate eigenvalues involved in the Hopf bifurcation has a singular limit as 
e — > [8 J . Guckenheimer [29] derives the following normal form for a singular Hopf bifurcation 

ex = y — x 2 , 

V = z-x, (18) 
z = —v — ax — by — cz, 

where x is a fast variable, (y, z) are slow variables and (u, a, b, c) are parameters. The key difference 
between ([15]) and (fT8j) is that we can find global equilibria q = q(u,a,b,c) for (fT8j) . They are 
determined by solving the equation 

- v = (a + c)x + bx 2 . (19) 

If v ~ then the equilibrium point is close to the folded singularity at the origin. The desingularized 
slow flow of (|18h can be calculated similar to the folded node case. It can be shown [151 50j that q is 
only important for the local dynamics near (0, 0, 0) if v is smaller than 0(e 1 / 2 ). The key differ cnce 
between MMOs that pass near a global equilibrium is that the SAOs can also be influenced by the 
stable and unstable manifolds W s (q) and W u (q). Detailed visualizations of the situation can be 
found in [15] [16] . Results for the unfolding of (fT8l) can be found in j29[ [32] . We are going to use 
the normal forms (|15f) and (|18p as "black-box" units for numerical simulation in Section [H 

B Error Analysis 

We briefly analyze the error of our approximation for the maps mx for K € {j, (a, +), b, /} for 
k = —10 and A G (^fsn, A n /)- The numerical integration of trajectories was carried out with a 
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Figure 13: .Affine and quadratic fit approximation error e of the different maps txik for K 6 
{j, (a, +), b, /} given in (|20|) . The horizontal axes are slices in A-parameter space with k = —10. 
The error is measured in three different norms L 1 (red), L 2 (green) and L°° (blue). The dashed 
curves for mj and rrif indicate the error for the lower branches and the solid curves the error for 
the upper branches of the maps. The domain for rrij and m a ^ + has been chosen as z G (2A — (4 + 
k) — 1, 2A — (4 + k) + 1) and the domain for and mj is the entire projection of the strong canard. 



standard stiff numerical integration method (ode 15s in MatLab [55]) with absolute error tolerance 
10~ 8 . The grid size h for the domain of the maps mx was always chosen so that h < 0.02. The 
main question we have to address is whether at a given fixed set of parameters (A, k) there exist 
affine and quadratic approximations as postulated in Section 2] Figure [13] shows the error of the 
fit to the postulated polynomial forms measured in three different norms 

<L l ) 
e(L 2 ) 
e(L°°) 

where m™ m indicates the map obtained from numerical integration and m 1 ^ 1 denotes the affine and 
quadratic fits. The integrals in (|20p have been evaluated from the discrete numerical integration 
data and the associated polynomials fits using a composite Simpson rule [65] which has error 0(h 5 ) 
as h — > 0. Figure [13] shows that the worst-case error for the proposed affine and quadratic maps 
due is at most on the order of 10 -2 over the entire range of parameters; the numerical integration 
error 10~ 8 and the numerical quadrature error h 5 < (0.02) 5 can be neglected here. Overall, the 
affine and quadratic approximations are certainly satisfactory to extract the basic MMO patterns. 

A natural question is to ask what happens to the perturbations of mx when e > 0. It is 
well-known from Fenichel theory that the error near normally hyperbolic segments of the critical 
manifold and in the fast subsystem is at most 0(e) as e — > 0. Near the fold points [68] it has 
been proven that the error is at most 0(e 1 ^ 3 ) as e — > 0. Therefore we have that rriK(z) + C(e 1//3 ) 
represents a flow map for < e <C 1. 

The numerical computations we present here can likely be made mathematically rigorous [36] 



= Kr um (*)-™£ (z)\dz, 

= {^{rn^ m {z)-rn^\z)fdz) l, \ (20) 
= snp ze[z0}Zl] \m^ m (z)-m^(z)\, 
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using interval arithmetic and tools such as IntLab [63]. The main reason for this conjecture is that 
rigorous numerical integration and quadrature are two standard situations in interval arithmetic 
|63j . However, carrying out this rigorous proof is beyond the scope and goals of this paper. 
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